
RNGkind("L'Ecuyer-CMRG")
set.seed(1234)


time.start <- proc.time()[3]

congs <- 101:113

load("all_data.Rda")
load("legDataTmp.Rda")

# Calculate the number of cores
no_cores <- min(detectCores() - 1, 10)

# define matrix inverse function
matrix_inverse <- function(X) chol2inv(chol(X))
  

#here we go

##set up the vote data
y <- spread(all_data[,c(2,3,6)], key = rollcall, value=vote)
rownames(y) <- y[,1]
y <- y[,-1]

##make speeches into a matrix
speeches <- spread(all_data[,c(2,3,5)], key = rollcall, value=leader_speech)
speeches <- speeches[,-1]
speeches <- as.matrix(speeches)
speeches[is.na(speeches)] <- 0

##Set things up
N <- nrow(y) #N legislators
J <- ncol(y) #J votes
C <- length(congs) #C congresses

##Initial values using Heckman-Snyder

y <- as.matrix(y)
lower <- dataset$legis.data$icpsrLegis[grep("PELOSI",rownames(dataset$votes))]
upper <- dataset$legis.data$icpsrLegis[grep("SENSENBR",rownames(dataset$votes))]
lower <- which(rownames(y)==lower)
upper <- which(rownames(y)==upper)

x.start <- tapply(2*all_data$pid-3,all_data$id,function(x)x[1])
al.start <- rep(0,J)
be.start <- rep(1,J)

x.current <- x.start
names(x.current) <- NULL
al.current <- al.start
be.current <- be.start
de.current <- matrix(0,N,C)

#we need to make a vector of time periods
times <- as.numeric(gsub("-.*","",colnames(y)))-min(as.numeric(gsub("-.*","",colnames(y))))+1
TimeMat <- matrix(0,J,C)
for (j in 1:nrow(TimeMat)) TimeMat[j,times[j]] <- 1

#hyperprior parameters
b0 <- rep(0,2)
B0 <- 25*diag(2)
t0 <- 0
T0 <- 1
d0 <- rep(0,length(unique(all_data$cong)))
D0 <- rep(3,length(unique(all_data$cong)))
dd0 <- matrix(c(t0,d0))
DD0 <- diag(c(T0,D0))


# #definte the speech-time matrix
legTimeMatrix <- vector("list", N)
for (i in 1:nrow(y)) legTimeMatrix[[i]] <- matrix(rep(speeches[i,],each=C),J,C,byrow=TRUE)*TimeMat

#define the update functions
ab.update <- function(j){
  solve(t(Xstar)%*%Xstar + solve(B0))%*%(t(Xstar)%*%Yhat[,j] + solve(B0)%*%b0)
}

ab.update_chol <- function(j){
  matrix_inverse(crossprod(Xstar) + matrix_inverse(B0))%*%(crossprod(Xstar,Yhat[,j]) + matrix_inverse(B0)%*%b0)
}


xd.update <- function(i){
  Xstar <- cbind(be.update,legTimeMatrix[[i]])
  solve(t(Xstar)%*%Xstar + solve(DD0))%*%(t(Xstar)%*%Yhat[i,] + solve(DD0)%*%dd0)
}

xd.update_chol <- function(i){
  Xstar <- cbind(be.update,legTimeMatrix[[i]])
  matrix_inverse(crossprod(Xstar) + matrix_inverse(DD0))%*%(crossprod(Xstar,Yhat[i,]) + matrix_inverse(DD0)%*%dd0)
}


##set tolerance holders
tol <- 0
count <- 0
deviance <- 100
llik <- NULL

##create a bunch of matrix holders
mu <- matrix(NA, N, J)
dmu <- matrix(NA, N, J)
pmu <- matrix(NA, N, J)


##start updating
while (tol<1-1e-6){
  
  
  ##create the NxJ matrix of legislator loyalty times speeches
  d_speech_mat <- as.matrix(de.current[,times]*speeches)
  
  ##calculate the matrix of means
  mu <- -matrix(rep(al.current,N),N,J,byrow=TRUE) +
    outer(x.current,be.current) +
    d_speech_mat
  
  #a minor correction for outliers
  mu[mu < -8] <- -8
  mu[mu > 8] <- 8
  
  #set the value of ystar for this iteration
  dmu <- dnorm(-mu)
  pmu <- pnorm(-mu)
  ystar <- (mu+dmu/(1-pmu))*(y==1) + (mu-dmu/(pmu))*(y==0)
  ystar[is.na(y)] <- mu[is.na(y)]
  
  ##set up variables and export for parallel estimation of roll call parameters
  Yhat <- ystar-d_speech_mat
  Xstar <- cbind(-1,x.current)
  
  ##estimate roll call parameters
  tmpRC <- mclapply(1:J, ab.update_chol, mc.cores=no_cores)
  tmpRC <- do.call(cbind,tmpRC)
  
  al.update <- tmpRC[1,]
  be.update <- tmpRC[2,]
  
  ##set up variables and export for parallel estimation of legislator parameters
  Yhat <- ystar + matrix(rep(al.current,N),N,J,byrow=TRUE)
  
  # Initiate cluster
  tmpLegis <- mclapply(1:N, xd.update_chol, mc.cores=no_cores)
  tmpLegis <- do.call(cbind,tmpLegis)
  
  x.update <- tmpLegis[1,]
  de.update <- t(tmpLegis[2:nrow(tmpLegis),])
  
  
  ##make sure orientation of x-axis is correct
  if (x.update[lower]>x.update[upper]){
    x.update <- -x.update
    be.update <- -be.update
  }
  
  
  ##calculate correlation-based model tolerance
  if (count>0) tol <- min(c(cor(x.update,x.current), cor(al.update,al.current),
                            cor(be.update,be.current), cor(c(de.update),c(de.current))))
  
  
  ##set new values as those just estimated
  x.current <- x.update
  al.current <- al.update
  be.current <- be.update
  de.current <- de.update
  
  count <- count + 1
  # print(count)
  
}


results <- list(x.current=x.current, be.current=be.current, de.current=de.current, al.current=al.current)
save(results, file=paste0(getwd(), "/EMestimates.Rda"))

cat("Replication file 1 done.","\n")


